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ABSTRACT 

The relationship between a decaying strong turbulence and kinetic instabilities in a slowly expand¬ 
ing plasma is investigated using two-dimensional (2-D) hybrid expanding box simulations. We impose 
an initial ambient magnetic held perpendicular to the simulation box, and we start with a spectrum 
of large-scale, linearly-polarized, random-phase Alfvenic fluctuations which have energy equipartition 
between kinetic and magnetic fluctuations and vanishing correlation between the two fields. A turbu¬ 
lent cascade rapidly develops, magnetic held fluctuations exhibit a power-law spectrum at large scales 
and a steeper spectrum at ion scales. The turbulent cascade leads to an overall anisotropic proton 
heating, protons are heated in the perpendicular direction, and, initially, also in the parallel direc¬ 
tion. The imposed expansion leads to generation of a large parallel proton temperature anisotropy 
which is at later stages partly reduced by turbulence. The turbulent heating is not sufficient to 
overcome the expansion-driven perpendicular cooling and the system eventually drives the oblique 
firehose instability in a form of localized nonlinear wave packets which efficiently reduce the parallel 
temperature anisotropy. This work demonstrates that kinetic instabilities may coexist with strong 
plasma turbulence even in a constrained 2-D regime. 


1. INTRODUCTION 

Turbulence in magnetized weakly collisional space and 
astrophysical plasmas is a ubiquitous nonlinear phe¬ 
nomenon that allows energy transfer from large to small 
scales, and, eventually, to plasma particles. Properties 
of plasma turbulence and its dynamics remain an open 


challen ging problem (Petrosyan et al. 
& Velli 1 2011). The solar wind constitu 


2010; [Matthaeus 
;es a natural lab¬ 


oratory tor plasma turbu lence (Bruno & Carbone 2013_ 
Alexandrova et ah||2013 ), since it offers the opportunity 


ot its detaifed diagnostics. Turbulence at large scales 
can be described by the magnetohydrodynamic (MHD) 
approximation, accounting for the dominant nonlinear 
coupling and for the presence of the ambient ma gnetic 
field t hat introduces a preferred direction (Boldyrev et al. 
2011). Around particle characteristic scales the plasma 
description has to be extended beyond MHD and, at 
these scales, a transfer of the cascading energy to par¬ 
ticles is expected. The solar wind turbulence indeed 
likely energizes particles: radial profiles of proton tem¬ 
peratures indicate an important heating which is often 
comparable to the estimated turbulent energy cascade 
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rate dM acBride et al. 2008 Cranmer et al. 2009 Hellinge 


et al. 2013 b This energization proceeds through coiii- 

sionless processes which may have a feedback on turbu¬ 
lence. In the solar wind the problem is further compli¬ 
cated by a radial expansion which induces an additional 
damping; turbulent fluctuations decrease due to the ex¬ 
pansion as well as due to the turbulent decay. T he expan¬ 


sion thus slows down the tur bulent cascade (cf., Grappin 
et al.|1 993 Do ng et al.|2014 )). Furthermore, the charac¬ 


teristic particle scales change with radial distance affect¬ 
ing possible particle energization mechanisms. 

Understanding of the complex nonlinear properties of 
plasma turbulence o n particle scales is facilitated via a 


nume rical approach (Servidio et al. 2015 Franci et al. 
2015a). Direct kinetic simulations of turbulence show 


that particles are indeed on average heated by the cas¬ 


cade (Parashar et al. 2009; Marlgr vskii k, Vasquez|[20TT 


Wu etal.||2013 ' |Franci et al.||2015a|), and, moreover, tur- 
bufence leads locally to compl ex anisotropic and nongy- 


rotropic distribution functions (Valentini et al.|2014 Ser¬ 


vidio et al. 2015). Furthermore, expansio n naturally gen¬ 


erate s particle temperature anisotropies (Matteini et al. 
2012). The anisotropic and nongyrotropic features may 


be a source of free energy for kinetic instabilities. In situ 
observations indicate existence of apparent bounds on 
the proton temperature anisotropies which are co nsiste nt 


with the oretical kinetic li near predictions (Hellinger et al. 
Hellinger & Travni'cek 2014). These line ar predic- 


2006 


tions have, however, many limite d assumptions (Matteini] 
et al.|2012 Isenberg et al.||2013 ); especially, they assume 
a homogeneous plasma which is at odds with the pres¬ 
ence of turbulent fluctuations. On the other hand, the 
observed bounds on the proton temperature anisotropy 
(and other plasma parameters) and enhanced magnetic 


fluctuations near these bounds (Wicks et al. 2013 La- 
combe et al.|2014~ ) indicate that tnese kinetic instabilities 


are active even m presence of turbulence. 
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2. SIMULATION RESULTS 

In this letter we directly test the relationship between 
proton kinetic instabilities and plasma turbulence in the 
solar wind using a hybrid expanding box model that al¬ 
lows to study self-consistently physical processes at ion 
scales. In the hybrid expanding box model a constant 
solar wind radial velocity v sw is assumed. The radial dis¬ 
tance R is then R = R 0 (l + t/t e o) where Ro is the initial 
position and t e0 = R 0 /v sw is the initial value of the char¬ 
acteristic expansion time t e = R/v sw = t e0 (l + t/t e o). 
Transverse scales (with respect to the radial direction) 
of a small portion of plasma, co-moving with the solar 
wind velocity, increase oc R. The expanding box uses 
these co-moving coordinates, approxim ating the spheri¬ 
cal coordinates by the Carte sian ones (Liewer et al.|200f 


Hellinger fc Travmcek|2005 ). The model uses the hybrk 


approximation where electrons are considered as a mass¬ 
less, charge neutralizing flui d and ions are described by a 
particle-in-cell model ([Matthews] 1994). Here we use the 
two-dimensional (2-D) version of the code, fields and mo¬ 
ments are defined on a 2-D x-y grid 2048 x 2048; periodic 
boundary conditions are assumed. The spatial resolution 
is A:r = Ay = 0.25d p o where d p o = vao/Q p o is the initial 
proton inertial length ( vao'- the initial Alfven velocity, 
flpo: the initial proton gyrofrequency). There are 1,024 
macroparticles per cell for protons which are advanced 
with a time step At = 0.05/U p o while the magnetic field 
is advanced with a smaller time step Atg = At/10. The 
initial ambient magnetic field is directed along the ra¬ 
dial, 0 direction, perpendicular to the simulation plane 
B o = (0,0, Bq) and we impose a continuous expansion 
in x and y directions. Due to the expansion the ambient 
density and the magnitude of the ambient magnetic field 
decrease as n oc B oc R~ 2 (the proton inertial length 
d p increases oc R, the ratio between the transverse sizes 
and d p remains constant; the proton gyrofrequency tt p 
decreases as a R~ 2 ). A small resistivity rj is used to 
avoid accumulation of cascading energy at grid scales; 
initially we set rj = 10~ 3 /.ioU^ o /D p o (yo being the mag¬ 
netic permittivity of vacuum) and y is assumed to be 
oc n. The simulation is initialized with an isotropic 2-D 
spectrum of modes with random phases, linear Alfven 
polarization (SB _L -Bo), and vanishing correlation be¬ 
tween magnetic and velocity fluctuation. These modes 
are in the range 0.02 < kd p < 0.2 and have a flat one¬ 
dimensional (1-D) power spectrum with rms fluctuations 
= 0.24Bo- For noninteracting zero-frequency Alfven 
waves the linear a pproximation predicts SB± oc R~ l 
(Dong et al. 2014]). Protons have initially the paral¬ 
lel^ proton beta j] = 0.8 and the parallel temperature 
anisotropy A p = T p ±/T p \\ = 0.5 as typical proton param- 
eters in the solar wind in the vi cinity of 1 AU (Hellinger 
et al.||2006| |Marsch et al.||2006| . Electrons are assumed 
to be isotropic and isothermal with /3 e = 0.5 at t = 0. 

The initial random fluctuations rapidly relax and a tur¬ 
bulent cascade develops. Figure [J shows the evolution of 
the 1-D power spectral density (PSD) Pg .l = -Pb_l(&) 
of the magnetic field B± perpendicular to B 0 . On large 
scales, the initial flat spectrum evolves to a power law. 
This large-scale power law remains clearly visible until 
t ~ 0.7t e o although its slope slowly varies in time, pass¬ 
ing from about -3/2 to -5/3 (these estimated slopes are, 
however, quite sensitive to the chosen range of wave vec¬ 


— t — 0.035t eO — t — 0.1t e o —t — 0.3t eO —t = 0.7t e0 —t — t e o 



Figure 1. (top) 1-D PSD Pbj l of the fluctuating magnetic field 
B l perpendicular to the ambient magnetic field and (bottom) Pb± 
compensated by /c 3 / 2 as functions of k at different times. The thin 
long dashed line shows the initial spectrum and the thin solid line 
shows a dependence oc A: 5/3 f or comparison. 


tors). The variation of large-scale slopes ( kd p < 1) is 
likely connected with the decay of the large scale fluctua¬ 
tions due to the cascade and the expansion as the inertial 
range is likely quite narrow. This problem is beyond the 
scope of the present letter and will be a subject of future 
work (note that a similar steepening is also obser ved in 
MHD expanding box simulations, cf., Dong et al.] |2014 ); 
this letter is mainly focused on ion scales. 

Around kd p ~ 1 there is a smooth transition in Pg± 
separating a the large-sc ale power law slo pe a nd a steeper 
slope at sub-ion scales (Franci et al. 2015b). The PSD 
amplitudes decay in time partly due to the expansion 
and partly to the turbulent damping. Note that there 
are some indications that the position of the transition 
shifts to smaller kd p with time/radial distance (compare 
the blue, green, orange, and red curves in Figure ll]); a 
similar trend is observed for the proton gyroradius since 
it increases only slightly faster than d p . At later times 
the fluctuating magnetic energy is enhanced at ion scales 
around ~ 0.4 4- 1 (compare the red and black curves 
in Figure IT]); this indicates that some electromagnetic 
fluctuations are generated at later times of the simula¬ 
tion. 

Figure[2]sunnnarizes the evolution of the simulated sys¬ 
tem which goes through three phases. During the first 
phase the system relaxes from the initial conditions and 
turbulence develops; the level of magnetic fluctuations 
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Figure 2. Time evolution of different quantities: (from top to 
bottom) the fluctuating magnetic field (solid) perpendicular \5B±\ 2 
and (dashed) parallel |5i3||| 2 with respect to Bq (the dotted line 
shows the linear prediction for the zero-frequency Alfven waves); 
the average squared parallel current (ji 2 ); the parallel T p y (solid 
line) and perpendicular Tp± (dashed line) proton temperatures 
(the || and _L directions are here with respect to the local magnetic 
field; the dotted lines denote the corresponding CGL predictions); 
(solid) the nonlinear eddy turnover time t n i at kd p = 1 (dotted) the 
expansion time t e , and (dashed) the linear time ti for the oblique 
firehose instability. 


increases at the expense of proton velocity fluctuations. 
The fluctuating magnetic field 5B±/B reaches the maxi¬ 
mum at about t ~ 0.008f e o. During this phjrse a parallel 
current j z is generated, (j?) normalized to B 2 /d 2 reaches 


a maximum at t ~ 0.035t e o indi cating a presence of a well 


developed turbulent cascade (Mininni & Pouquet 2009 


Valentini et al. 2014). After that the system is dominatec 


by a decaying turbulence, the fluctuating magnetic field 
initially decreases faster than B till about 0.3f e o- 
During the second phase protons are heated. For neg¬ 
ligible heat fluxes, collisions, and fluctuations, o ne ex¬ 
pects the double adiabatic behavior or CGL (|Chew et~ah 


1956 Matteini et al.|2012|: the parallel and perpendicu¬ 


lar temperatures (with respect to the magnetic field) are 
expected to follow T p ± oc B and T p y = const., respec¬ 
tively. T p l decreases slower than B during the whole 
simulation, protons are heated in the perpendicular di¬ 
rection while in the parallel direction the heating lasts till 


about t ~ 0.25f e o whereas afterwards protons are cooled. 
The parallel and perpendicular heating rates could be 
estimated as (cf., Verscharen et al.|[2015|: 


n 3 fcs d (T\\B 


B 2 d t 


d (T 

and Q± = nk B B — ( — 
at \ B 


Q ii = 

J_J KXO \ 10 / \ JLJ J 

( 1 ) 

A more detailed analysis indicates that between t = 
0 . 1 t e o and t = 0.7f e o the parallel heating rate Q || 
smoothly varies from about 0 . 2 Q e and — 0 . 2 Q e , whereas 
Q± is about constant ~ 0.2Q e ; here Q e = nkBT/t e . In 
total, protons are heated till t ~ 0.7£ e o and the heat¬ 
ing reappears near the end of the simulation t > 0.95f e o. 
Note that the perpendicular heating rate is a nonnegli- 
gible fraction of that observed in th e solar wind where 
<3_l ~ 0.6Q e (Hellinger et al. 2013); however, the pro¬ 
ton heating in 2D hybrid simulations is typic ally quite 
sensitive to the used electron equation of state (Parashar 
et al.|2014 ) and also to the used resistivity a nd the num¬ 
ber of particles per cell ( Franci et al.| 2015a). The turbu¬ 
lent heating is, however, not sufficient to overcome the 
expan sion-driven perpendicu lar cooling as in the solar 
wind (Matteini et al. 2007). During the third phase, 
t > 0.7t e o, there is an enhancement of the parallel cool¬ 
ing and perpendicular heating which cannot be ascribed 
to the effect of the turbulent activity. For a large par¬ 
allel proton temperature anisotropy a firehose instability 
is expected. The presence of such an instability is sup¬ 
ported by the fact that the fluctuating magnetic field 
increases (with respect to the linear prediction) suggest¬ 
ing a generation of fluctuating magnetic energy at the 
expense of protons. To analyze the role of different 
proces ses in the system we e stimate their characteristic 
times ( |Matthaeus et al.|2014 |. The bottom panel of Fig¬ 
ure [ 2 ] compares the turbulent nonlinear eddy turnover 
time t„,i = fc~ 3 / 2 (P R | (fc)//inutn )~ 1 / 2 at kd p = 1 (cf. 


Matthaeus et al. 2014), the expansion tim e t P , and the 
linear time ti of th e oblique firehose (Hellinger & Mat-| 
sumotopOOO 2001 1 estimated as ti = l/ 7 m, where 7 m is 
the maximum growth rate calculated from the average 
plasma properties in the box assuming bi -Maxwellian 
proto n velocity distribution functions (Hellinger et al. 
2006). The expansion time t e is much longer than t n i at 
kdp = 1 (as well as at the injection scales). The expand¬ 
ing system becomes theoretically unstable with respect 
to the oblique firehose around t ~ 0.47f e o but clear signa¬ 
tures of a fast proton isotropization and of a generation of 
enhanced magnetic fluctuations appear later t > 0.7t e o. 
This is about the time when the linear time becomes 
comparable to the nonlinear time at ion scales. After 
that, tiklp slightly increases as a result of a saturation of 
the firehose instability whereas t n ifl p at kd p = 1 is about 
constant (note that kl p decreses as R~ 2 ). This may indi¬ 
cate that the instability has to be fast enough to compete 
with turbulence; however, the 2-D system has strong geo¬ 
metrical constraints. Also the stability is governed by the 
local plasma properties. Figure [3] shows the evolution of 
the system in the plane (/3 p ||, Ap). During the evolution, 
a large spread of local values in the 2-D space (0 P \\, A p ) 
develops. Between t — 0.1t e o and t — 0.65t e o the aver¬ 
age quantities evolve in time following (A p ) oc (/3 p ||) -0 ' 86 - 
This anticorrelation is qualitatively similar to in situ He- 
lios observations between 0.3 and 1 AU (Matteini et al. 
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2007). During the third stage, when the strong parallel 


temperature anisotropy is reduced, both local and aver¬ 
age values of fi p \\ and A p appear to be bounded by the 
linear ma rginal stabi lit y condition s of the oblique fire¬ 
hose (Hellinger & Travnfcek 2008), although relatively 


large theoretical growth rates ~ O.lfL, are expected. 



Figure 3. Simulated data distribution in the plane (P p \\,Ap) at 
different times. The empty circles give the initial condition whereas 
the solid circles denote the average values. The solid lines show 
the evolution of the average values. The dashed contours show the 
maximum growth rate 7 m (in units of Q p ) of the oblique firehose 
instability as a function of j3 p \\ and A p . The dotted contours display 
the corresponding angle of propagation of the most unstable mode. 


to the magnetic field near threshold, further away from 
threshold the unstable modes become more oblique (see 
Figure [3]) and these oblique angles are locally available 
between vortices where we observe the enhanced level of 
magnetic fluctuations. These geometrical factors may be 
responsible for the late appearance of the instability (but 
constraints on the instability time scales imposed by the 
turbulent non-linearities are likely also important). The 
localized wave packets are Alvenic, have wavelengths of 
the order of 10d p and propagate with a phase velocity 
about 0.1 va in agreement with expectation for the non¬ 
linear phase of the oblique firehose. Furthermore, the 
parallel temperature anisotropy is strongly reduced in 
their vicinity. These Alvenic wave packets are responsi¬ 
ble for the enhanced level of the magnetic PSD at ion 
scales seen in Figure [l] 

For a linear instability it is expected that the mag¬ 
netic fluctuations increase exponentially in time during 
its initial phase (except when th e growth time is compa - 
rable to_the expansion time, cf., Tenerani & Vclli 2013). 
Figure [2] however shows that the overall magnetic fiuc- 
tuations 5B±/B (with respect to the linear prediction) 
and 6B\\/B increase rather slowly (secularly) in time for 
t > 0.7 t e . This behavior is expected for a long ti me evolu- 
tion in a fo rced system a fter the saturatio n (cf., Matteini] 


et al.||2006 Rosin et al.||2011 Kunz et al.||2014 ). An ad¬ 

ditional analysis indicates that the expected exponential 
growth is indeed seen in the simulation but only locally 
both in space and time. This exponential growth is ob¬ 
scured by the turbulent fluctuations and, furthermore, it 
is blurred out due to the averaging over the simulation 
box in the global view of Figured] 

On the microscopic level the firehose activity leads to 
an efficient scattering from parallel to perpendicular di¬ 
rection of protons in the velocity space. Figure [5] shows 
the evolution of the proton velocity distribution function 
/ = f{v\\,v±) averaged over the simulation box. While 
turbulenc e leads locally to a complex proton distribu tion 
functions ( |Valentini et al.|2014 Servidio et al.|20l5 cf.,) 
the average proton distribution function during the first 
two phases remain relatively close to a bi-Maxwellian 
shape (Figure |5j top panel). During the third phase 
there appear clear signatures of the cyclotron diffusion 
(for protons with f || > va) as expected for the oblique 
firehose instability (Hellinger & Travnfcek 2008). 


The evolution in the real space is shown in Figure [4] 
which shows the magnitude of the perpendicular fluc¬ 
tuating magnetic field 5B j_ and the proton temperature 
anisotropy A p at different times (see also the movie which 
combines the evolution in Figures [3] and [4]). The modes 
with initially random phases rapidly form vortices and 
current sheets. Despite the overall turbulent heating a 
strong temperature anisotropy T p ± < T p y develops ow¬ 
ing to the expansion and a firehose-like activity devel¬ 
ops in a form of localized waves/filaments with enhanced 
SB±. These fluctuations appear in regions between vor¬ 
tices where B±/B is enhanced, i.e., in places where the 
angle between the simulation_ plane and the local mag¬ 
netic field 9b arccos B±/B) is less oblique (reaching 
~ 60° and below). This is in agreement with the the¬ 
oretical expectations, while the oblique firehose is un¬ 
stable for moderately oblique wave vectors with respect 


3. DISCUSSION 

Using 2-D hybrid simulations we investigated the evo¬ 
lution of turbulence in a slowly expanding plasma. The 
numerical model shows that the turbulent heating is not 
sufficient to overcome the expansion driven cooling and 
that the oblique firehose becomes active for a sufficiently 
large parallel proton temperature anisotropy and for suf¬ 
ficiently oblique angles of propagation. 

While the modeled expansion is about ten times faster 
than in the solar wind, the ratio between the expansion 
and the nonlinear eddy turnover time scales is quite re¬ 
alistic: t e /tni ~ 1000 at kd p = 1 for t > 0.7t e which is 
about four times smaller than that of the solar wind with 


simila r plasma parameters at 1 AU (Matthaeus et al. 
2014|. Note also that a similar evolution is observed tor 


many different plasma and expansion parameters. 

In the present case, both turbulence and the 2-D ge- 
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Figure 4. Color scale plots of (left) SBj_ and (right) A p as func¬ 
tions of x and y for (top) t = 0.1t e o, (middle) t = 0.7f e o, and 
(bottom) t = t e o- The solid lines show selected (projected) mag¬ 
netic field lines. Only a quarter of the simulation box is shown. 


ometry constraints strongly affect the firehose instability 
and there are indications that firehose has an influence 
on turbulence (the m ixed third-order structure functions 
(Verdini et al. 2015) are enhanced due to the firehose 
activity suggesting a stronger cascade rate). The prob¬ 
lem of the interaction between turbulence and kinetic in¬ 
stabilities requires further work; three-dimensional sim¬ 
ulations are needed to investigate the interplay between 
turbulence and instabilities as usually the most unsta¬ 
ble modes are parallel or moderately oblique with re¬ 
spect to the ambient magnetic field; in the present case 


the p arallel firehose (Gary et al. 


20061 would be the dominant instability but the 2-D con- 


1998 Matteini et al. 


straints strongly inhibit it. On the other hand, numer¬ 
ical simulations indicate that the oblique firehose plays 
an important role in constraining the proton tempera¬ 
ture anisotropy in the expanding solar wind even in the 


case when the p arallel firehose is dominant (Hellinger & 
Travnfcek 2008). Nevertheless, the present work for the 
first time clearly demonstrates that kinetic instabilities 
may coexist with strong plasma turbulence and bound 
the plasma parameter space. 


The authors wish to acknowledge valuable discussions 
with Marco Velli. The research leading to these results 
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Figure 5. Average proton velocity distribution function / as a 
function of parallel and perpendicular velocities un and v± (with 
respect to the local magnetic field) for for (top) t = 0.1t e o, (middle) 
t = 0.7f e o, and (bottom) t = t e q. 
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